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Abstract 

The  immersed  finite  element  method  based  on  a  uniform  Cartesian 
mesh  has  been  developed  for  elasticity  equations  with  discontinuous  phys¬ 
ical  parameters  across  an  interface  in  this  paper.  The  interface  does  not 
have  to  be  aligned  with  the  mesh.  The  main  idea  is  to  modify  the  basis 
function  over  those  triangles  in  which  the  interface  cuts  through  so  that 
the  natural  interface  conditions  are  satisfied.  The  standard  linear  basis 
functions  are  used  for  other  triangles.  A  level  set  function  whose  zero  level 
set  represent  the  interface  is  used.  Numerical  examples  are  also  presented. 
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1  Introduction 


In  this  paper,  we  develop  a  finite  element  method  using  a  uniform  Cartesian 
mesh  for  a  two-phase  elasticity  system  of  equations  of  the  form: 
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v  •  cr  +  F  =  0  in  n,  (1.1) 

[  u  ]  =  0  on  T,  (1.2) 

[  crn  }  =  0  on  T,  (1.3) 

u  =  u0  on  Oil,  (1.4) 

where  cr  is  the  stress  tensor,  F  =  (Fi,F2)t  :  O  — >  K2  is  the  body  force, 
u  =  (u\,U2)t  is  the  displacement  field,  T  is  a  smooth  interface  that  divides 
the  domain  O  into  two  parts  0+  and  fi~,  n  =  {n\,n2)T  is  the  unit  normal 
vector  of  the  interface  T,  pointing  from  the  —  phase  to  the  +  phase,  and  u0  is  a 
given  vector-valued  function  that  represents  the  displacement  on  the  boundary 
dil.  Across  the  interface  T,  the  physical  parameters  such  as  Young’s  modulus 
and  Lame  constants  have  a  finite  discontinuity.  For  a  function  v,  We  use  [i>]  = 
v+  —  v~  to  denote  the  jump  of  v  across  the  interface  T.  The  jump  condition 
(1.2)  means  that  u  is  continuous  across  the  interface. 

Multi-phase  elasticity  problems  often  arise  in  materials  science,  see  for  ex¬ 
ample,  [12,  15].  We  refer  the  reader  to  [1,  6,  9,  17]  for  some  applications  and 
the  references  therein.  However,  solving  such  elasticity  system  is  often  difficult 
due  to  the  arbitrary  interface  and  the  discontinuities  in  the  coefficients  and  the 
gradient  of  the  solution. 

There  exist  several  numerical  methods  for  solving  general  elasticity  problems 
that  do  not  involve  interfaces.  Among  them,  the  finite  element  method  and  the 
boundary  integral  or  boundary  element  method  appear  to  be  very  successful,  cf. 
e.g.,  [2,  13,  18]  and  the  references  therein.  However,  in  treating  moving  interface 
problems,  the  use  of  fixed  Cartesian  grids  often  shows  advantages  in  practical 
computations  [7].  It  is  therefore  desirable  to  develop  new,  efficient  methods 
based  on  fixed  Cartesian  grids  for  the  elasticity  system  with  an  interface. 

In  [17],  a  second  order  immersed  interface  method  was  developed  for  the 
elasticity  system  with  an  interface  based  on  Cartesian  grid.  However,  due  to 
lack  of  the  maximum  principle,  the  stability  of  the  method  is  still  under  investi¬ 
gation,  and  the  linear  solver  for  the  finite  difference  equations  may  not  converge 
satisfactorily.  The  goal  of  this  paper  is  to  develop  an  immersed  finite  element 
(IFE)  for  solving  the  two-phase  elasticity  system.  The  idea  using  a  Cartesian 
mesh  to  solve  a  single  elliptic  interface  problems  can  be  found  in  [10,  11]. 

This  paper  is  organized  as  follows:  In  Section  2,  we  derive  the  weak  form  of 
the  elasticity  system  and  the  outline  of  the  finite  element  method.  In  Section  3, 
we  discuss  how  to  constructs  the  basis  functions  for  non-interface  and  interface 
elements.  In  Section  4,  we  explain  how  to  use  a  the  level  set  function  to  represent 
the  interface.  Several  numerical  examples  are  given  in  Section  5. 
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2  Derivation  of  the  weak  form  the  elasticity  sys¬ 
tem 


In  this  section,  we  derive  the  weak  formulation  of  the  elasticity  system.  We  use 
more  general  boundary  conditions  (2.1)  and  (2.2)  which  replace  (1.4) 

crn  =  t  ,  on  (2-1) 

u  =  u o  ,  on  (2-2) 


where  <9S2  =  dQi  U  d^-  For  convenience,  we  rewrite  the  strain  and  stress  as 
vector  forms 


r  _ 
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du 
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£x 

dx 
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If  we  introduce  operator  A  as  follows 
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then  the  strain-displacement  and  stress-strain  relations  can  be  rewritten  as 
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(2.3) 


O"., 

&X 

(Ty 

=  [D].e=[D]. 

Sy 

-  Txy  - 
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(2.4) 


where  [  D  ]  is  the  elasticity  matrix  (or  constitutive  stress-strain  matrix), 


[D] 


A  +  2  n  X  0 

A  A  +  2^  0 

0  0  n 


(2.5) 


in  which  A  and  /i  are  Lame  constants.  If  we  let  E  be  the  Young’s  modulus,  and 
v  be  the  Poisson’s  ratio,  then  we  have 


A  = 


Ev 
1  —  v2 


E 

^  2(1+  vy 


(plane  stress), 


A  = 


Ev 

(1  —  2  v)(l  +  v) 


(plane  strain). 
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The  equilibrium,  constitutive  and  strain-displacement  equations  then  become 


ATcr  =  —F, 

(2.6) 

a  =  [D]e, 

(2.7) 

£  =  A  u  . 

(2.8) 

Eliminating  er  and  e  gives  the  “displacement”  formulations 

At[  D]Au  =  -F. 

(2.9) 

Equation  (2.9)  will  be  used  to  derive  the  stiffness  matrix. 


2.1  The  variational  form 

Let  us  consider  the  potential  energy  II  of  an  elastic  body.  II  is  defined  as  the 
sum  of  the  total  strain  energy  ( U )  and  the  work  potential  (WP) 

II  =  Strain  energy  +  Work  Potential. 


For  linear  elasticity  materials,  the  strain  energy  per  unit  volume  in  the  body  is 
\  crT  e  .  For  a  elastic  body,  the  total  energy  U  is  given  by 

U  =  -  f  crT  e  dtt  =  -  I  £T[  D  ]  £  dfl. 

2  Jn  2  Jn 

The  work  potential  is  given  by 

wp  =  -  f  utf  dn-  uTt  ds. 

J  J  Qoi 

The  total  potential  for  the  general  elastic  body  is 

[  uTt  dS.  (2.10) 

fioi 

By  the  principle  of  minimum  potential  energy  (cf.  [3,  4,  5,  8,  14]),  we  obtain  the 
following  “variational  form”  or  “weak  form”  for  two-dimensional  stress  analysis 
(  m  is  an  arbitrary  displacement  field): 

[  <jt  £  (  u  )  dV  =  I  uTf  d.n  +  [  uTt  dS.  (2-11) 

J 17  J  £1  J  Qoi 


D  }  £  dO  u1  F  dil  — 
Jn 


2.2  The  basic  finite  element  equations 

We  assume  the  domain  O  is  a  rectangle,  but  the  interface  T  can  be  arbitrary. 
We  use  a  uniform  triangulation  regardless  of  the  interface  T.  Therefore,  the 
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interface  generally  is  not  aligned  with  the  edges  of  the  triangulation.  We  need 
to  find  the  basis  functions  for  each  right  triangle  of  the  partition.  For  each 
triangle,  let  {  u  *}  be  the  vector  of  the  nodal  displacements: 

{  u  *}  =  {«i  Vi  u2  v2  u3  1>3}. 


The  displacements  at  a  point  inside  the  triangle  u  can  be  determined  in  terms 
of  the  nodal  displacements  {  u  *}  and  the  basis,  or  shape,  functions  N  : 

u  =  [  N  ]{  u*  }. 


Strains  and  stresses  can  also  be  determined  at  nodal  displacements: 

e  =  Au  =  A  [  N  ]{  u*  }  =  [  B  ]{  u*  }, 

<r  =  [£>]  e  =  [D][  B  ]{  «*  }, 

where  [  B  ]  =  A  [  N  ]  is  called  the  displacement  differentiation  matrix,  which 
can  be  obtained  by  differentiating  displacements  expressed  through  shape  func¬ 
tions  and  nodal  displacements.  So  we  can  obtain  the  following  equilibrium 
equations  for  a  finite  element: 

[  [B]t[D][B  u  }  =  f  [  N  }T  F  dxdy  +  [  [  N  }T  t  dS. 

JQ  Jq,  J  dQi 

On  each  element,  we  have 


/ 


[B]T[D][B}dn{u}=  [  N  }T  F  dxdy  + 


/  A(e) 


Ian  inA<e) 


[  N  }T  t  dS. 


We  need  to  know  the  basis  function  N  for  each  element. 


3  Constructing  the  basis  functions  at  interface 
triangles 

The  triangles  in  our  partition  are  classified  into  two  categories:  the  interface 
triangles  if  the  interface  divides  the  triangle  into  two  subsets,  and  non-interface 
triangles  otherwise.  For  non-interface  triangles,  we  use  standard  linear  basis 
functions.  However,  for  interface  triangles,  we  use  an  undetermined  coefficients 
method  to  determine  the  basis  functions  by  enforcing  the  natural  interface  con¬ 
ditions.  In  the  following  discussions,  we  discuss  how  to  construct  such  piecewise 
linear  basis  functions  for  interface  triangular  elements. 

For  a  typical  interface  triangle  A  ABC,  let  B  =  (#1,2/1),  C  =  (#2,2/2),  A  = 
(#3,2/3)  be  its  vertices.  Let  D  =  (#d,2/d)  and  E  =  (xe,ye)  be  the  intersections 
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A{x3,  ») 


Figure  1:  A  typical  interface  triangle  A  ABC. 


of  the  interface  and  the  edges  of  the  triangle,  see  the  sketch  in  Figure  1  for  an 
illustration.  For  convenience,  we  denote  T+  the  region  above  the  interface,  and 
T~  the  region  below  it. 

Once  the  values  u,  v  at  vertices  A,  B ,  C  of  the  element  T  are  specified,  we 
construct  the  following  piecewise  linear  function 


u(x,y) 
v(x,  y) 


u+(x,y)  =  ai  +  a2x  +  a3y,  if  (i,j)er+, 
u~{x,y)  =  bi  +  b-2x  +  b3y,  ii{x,y)  &T~\ 

v+(x,y)  =  ci  +  c2x  +  c3y,  if  (x,y)eT+, 
v~{x,y)  =  di  +  d2x  +  d3y,  if  (x,  y)  £  T~ , 


(3.1) 

(3.2) 


where  a*’s,  bi  s,  Cj’s,  d,;’s,  (i  =  1,2,3)  are  undetermined  coefficients.  The  values 
of  u,  v  at  the  vertices  A,  B1  C  are  u+(A)  =  u\,  u~(B)  =  u2,  u~(C )  =  u3, 
u+(A)  =  Vi,  v~(B)  =  v2,  v~(C)  =  v3.  Therefore  we  have 


-4(:E3,  y3)  :  | 

a3  +  a2x3  +  a3y3  —  u3 , 

,  Cl  +  c2x3  +  c3y3  =  v3\ 

(3.3) 

B(xuyi)  ■  j 

f  bi  +  b2Xi  +  b3yi  =  Mi, 
l  di  +  d2a"i  +  d3yi  =  v3\ 

(3.4) 

C(x2,y2)  :  j 

\  h+b2x2  +  b3y2  =  u2 , 

L  d  1  +  d2x 2  +  d3y2  =  v2\ 

(3.5) 

D(xd,  yd)  ■  < 

f  ai  +  a2xd  +  a3yd  -  bi  -  b2xd  -  b3yd  = 
t  Cl  +  c2xd  +  c3yd  -  di-  d2xd  -  d3yd  = 

cT  0 

(3.6) 

E(xe,  ye)  :  j 

01  +  a2xe  +  a3ye  -  bi  -  b2xe  -  b3ye  = 

[  Cl  +  c2xe  +  c3ye  -  di  -  d2xe  -  d3ye  = 

0  0 

(3.7) 
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By  the  interface  conditions, 


[  an  ]  =  0  on  the  interface  T, 


where  n 


n  i 
n2 


dicular  to  DE,  and 


is  the  unit  normal  vector  of  the  interface,  which  is  perpen- 


an 


(Tuni  +  o-i2n2 
a  i2?7-i  +  cr22n2 


In  the  component  form,  we  have 


du  dv  ( du  dv\ 

(A  +  2„) -n1  +  X-n1  +  ^-  +  -J  n2 ^ 

f  du  dv\  du  dv 

f‘(^  +  &J"1  +  A*:’‘2  +  (A  +  2'‘)%"2 

where  (assuming  the  plane  deformation) 

E 


=  0, 


=  o, 


d  = 


2(1  +  v)  ’ 


and 


A  = 


Ev 


(1  +  i/)(l  —  2v) 


(3.8) 

(3.9) 


If  we  use  [  G  ]  to  represent  the  coefficient  matrix  for  the  equations  (3.3)  ~  (3.9), 
and  use  X  to  represent  the  vector  formed  by  a\,  a2,  03,  61,  fe2,  63,  ci,  c2,  C3,  d\, 
d2,  c?3,  then  (3.3)  ~  (3.9)  can  be  rewritten  as  a  matrix-vector  form  as  follows: 


[G]X  =[C]{u*}, 


where 


0  0  0  0  1  0 
0  0  0  0  0  1 
1  0  0  0  0  0 
0  1  0  0  0  0 
0  0  1  0  0  0 
0  0  0  1  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 


{«*} 


Solving  3.10  gives 

X  =  [  G  ]-1[  C  ]{  m  *}. 


u  1 

Vi 

U2 

V2 

U3 

V3 


(3.10) 


(3.11) 
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In  other  words,  from  (3.3)  ' 
ai,  a2,  a3,  bi,  b2,  b3,  a,  c2 
{  u  *}.  Therefore,  u+,  u~ , 

(3.7),  and  the  interface  conditions,  we  can  express 
,  c3,  di,  d2,  d3  in  terms  of  Mi,  Vi,  u2,  v2,  u3,  v3,  or 
m+,  t>+,  can  be  expressed  as 

u+  = 

(  a1)  +  x~a2  +  y  a3) 

•{«*}, 

(3.12) 

u~  = 

(  bi  +  x  b2  +  y  b3  ) 

•{«*}. 

(3.13) 

v+  = 

[~Ci  +x~c^  +  y~c, i4)  • 

{«*}, 

(3.14) 

v~  = 

^  d±  +  x  d2  +  y  d3  ^ 

•{«*}, 

(3.15) 

where  a  (  =  (an,  ai2,  ai3,  «i4,  a^,  aie)  is  the  first 
in  (3.11),  and  so  forth  for  a  2,  •  •  • ,  d3.  If  we  set 

row  of  matrix  [  G  ]“ 

-Me] 

u 

'  Nn 

Nu 

N 13 

N 14 

N 15 

N lg 

V 

Nil 

n22 

n23 

A24 

n23 

N2Q 

Ml 

Ml 

U2 

V2 

U3 

V3 


def 


[N  ]{«*}, 


where  Ntj 's  ( i  =  1, 2,  j  =  1,  2,  •  •  •  ,  6)  are  piecewise  linear  functions,  then  we 
have 


ai  +  x  a2  +  y  a3 

~c±  +  x~c?  +  y~c^ 

_  _ ^  >  >  _ 

bi  +  x  b2  +1/^3 

di  +  x  d2  +y  d3 

+  side, 


—  side; 


(3.16) 


Figures  2  and  3  show  the  mesh  and  contour  plots  of  a  pair  of  the  basis 
functions  in  an  interface  triangles  and  in  its  entire  support. 

We  then  can  compute  the  differentiation  strain  matrix  [  B  }  for  interface 
elements: 


[B}+  =  A  [AT] 


a3 


c2 


[B}~ 


63  +  d2 


(3.17) 


A  [AT] 


(3.18) 


Shape  function  for  u 


Shape  function  for  v 


Figure  2:  The  mesh  and  contour  plots  of  a  pair  of  local  basis  functions  ( u ,  v) 
over  an  interface  triangle.  The  parameters  are  A+  =  80,  A-  =  160,  v+  =  0.35, 
and  v~  =  0.15. 


A  global  basis  function  Contour  plot 


Figure  3:  The  mesh  and  contour  plots  of  a  global  basis  function  u  on  its  support. 
The  parameters  are  A+  =  40,  A~  =  90,  v+  =  0.35,  and  =  0.15. 
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Therefore,  the  local  stiffness  matrix  is 


[  km  ] 


I  [B]T[D][B]dxdy 

M*cl 

f  [  B +]T[D +][  B +}dxdy 

J  trianale 


(3.19) 


+  f  [B  ~]T[  D  "][  B  -}dxdy 

J other 

[  B  +]T[  D  +][  B  +]  ■  (the  area  of  the  triangle) 

+  [  B  ~]T[  D  ”][  B  “]  •  (the  area  of  the  other  part), 


which  is  also  a  constant  and  symmetric  matrix. 


4  The  interface  and  the  level  set  method 


We  use  the  zero  level  set  of  a  Lipschitz  continuous  function  <f>(x,  y)  to  represent 
the  interface  T.  Usually  <t>{x,  y)  is  chosen  as  the  signed  distance  function.  With 
the  representation  of  a  level  set  function,  it  is  easy  to  compute  the  geometric 
information  that  is  needed  to  construct  the  basis  functions. 


Theorem  4.1.  Assume  the  coordinates  of  A,  B,  C,  D,  E  ( see  Figure  4)  are 
A(xi,yj+ 1),  B(xi,yj),  C(xi+i,yj),  D(xD,yD),  E(xe,Ve)-  Then  we  have 


xD 

l Id 


xE 


Ue 


Xi, 
Vj  + 

Vi+ 1 


fajAy 

'■  ,j  ~  4>i,j+ 1 


Ax(fii,j+ 1 


0((  Ay)2) 


0((  Ay)2), 


X-i+l  + 


-A  xcj>x(E)  +  A  y(j>y(E) 


-A  xcj>x(E)  +  A  yf>y(E) 


0(Aa;Ay) 

+  O(AxAy), 


Ay 

Vi  +  ^(*i+i  -  xE), 


where  1)?  ^('U+i, yj)>  and 


(4.1) 

(4.2) 


(4.3) 


(4.4) 


ME)  = 

M  Mij+0{AxAy)j 

(4.5) 

ME)  = 

M+i  -  M  +  o(AxAy). 

(4.6) 

Vj+ 1  -  Vj 
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The  above  theorem  gives  all  the  information  that  is  needed  to  construct  the 
basis  function  on  an  interface  triangle. 


Figure  4:  The  geometric  information  of  an  interface  triangles  and  the  level  set 
function  representation. 

The  proof  of  the  theorem  is  straightforward  using  the  Taylor  expansion  at 
the  intersections  D  and  E. 


5  Numerical  examples 

In  this  section  we  present  three  numerical  examples.  Without  loss  of  generality, 
the  domain  Q  is  a  rectangle  [— 1, 1]  x  [— 1, 1],  and  the  interface  T  is  x2  +  y1  = 

In  Example  5.1  andEXAMPLE  5.2,  we  consider  the  system  with  zero  body  forces 
and  homogeneous  jumps  conditions  for  displacements  and  traction  across  the 
interface.  So  the  system  is  as  follows, 


Vo-  = 

0, 

in  fl+  U  Q 

[«]  = 

0, 

on  r, 

<771  ]  = 

0, 

on  r, 

U  = 

uo  , 

on  dfl, 

where  Uq  =  [uo(x,y),  Vo(x,y)]T.  Since  the  analytic  solution  is  not  available, 
we  compare  the  computed  solutions  obtained  from  the  immersed  finite  element 
method  with  the  solutions  obtained  from  the  immersed  finite  difference  method 
developed  in  [16,  17]  for  example  5.1  to  validate  the  IFEM  method  for  the  elas¬ 
ticity  systems  with  interfaces.  The  immersed  finite  difference  method  developed 
in  [16,  17]  have  been  tested  against  the  exact  solution.  Example  5.2  is  similar 
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to  Example  5.1  except  very  large  contrast  (large  jump  ratio)  of  the  physical 
parameters  in  which  the  finite  difference  method  converges  slower  than  the  finite 
element  element  approach.  In  the  last  example  Example  5.3,  we  show  a  test 
result  with  non-zero  source  terms. 


Example  5.1.  In  this  example,  the  differential  equations,  interface  and  bound¬ 
ary  conditions  are  exactly  the  same  as  that  in  Example  2.5.3  in  [16]. 
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u0{x,y)  =  (x2  +  y2)2  +  —  In  (2  x1  +  y2) 

Do (x,y)  =  In  (l  +  x2  +  3  y2)  +  sin  (xy)  —  4x2  —  4  y2  +  1. 

The  parameters  are 


v~  =  0.22,  /Jr  =  200,  E~  =  488, 
v+  =  0.11,  y+  =  60,  E+  =  133.20. 


Figure  5  and  Figure  6  are  the  mesh  and  contour  plots  of  the  solution  (u,v) 
obtained  by  the  IFEM  method  described  in  this  paper  using  a  56  by  56  grid, 
which  generates  the  system  of  size  6272  x  6272.  From  Figure  5  and  Figure  6,  we 
can  see  they  are  almost  identical  to  the  figures  by  the  finite  difference  method 
in  Example  2.5.3  in  [16].  The  finite  difference  method  is  a  little  than  the  finite 
element  method. 


Figure  5:  The  mesh  and  contour  plots  of  the  computed  solution  u  of  Exam¬ 
ple  5.1  obtained  by  the  immersed  finite  element  method  developed  in  this  paper. 

Example  5.2.  In  this  example,  we  use  the  same  differential  equations,  the  jump 
conditions  as  that  in  Example  5.1.  But  we  use  different  boundary  conditions 
and  elasticity  constants  as  follows: 


=  0.02, 

=  0.4902, 

E~  = 

1, 

v+  =  0.49, 

=  167.7852, 

E+  = 

500, 

uo(x,  y) 

=  cos 

+  1 

f)+y  +  1  + 

sin  ((x 

+  1)  (tf+1)) 

vo{x,y) 

=  (x  + 1) 

(?/  +  1)  +  sin  (2  (x  +  1)  (y  +  1)) 
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Figure  6:  The  mesh  and  contour  plots  of  the  computed  solution  v  of  Exam¬ 
ple  5.1  obtained  by  the  immersed  finite  element  method  developed  in  this  paper. 


We  have  large  jump  ratio  in  the  parameters.  Figure  (7)  and  Figure  (8)  are  the 
mesh  and  contour  plots  for  the  displacements. 

In  this  example  the  finite  difference  method  takes  much  longer  time  to  con¬ 
verge  than  the  finite  element  method  does.  This  is  because  the  linear  system  of 
equations  obtained  from  the  finite  element  method  is  symmetric  and  has  bet¬ 
ter  condition  number  compared  with  that  obtained  from  the  finite  difference 
approach. 


Figure  7:  The  mesh  and  contour  plots  of  the  computed  solution  u  of  Exam¬ 
ple  5.2  obtained  by  the  immersed  finite  element  method  developed  in  this  paper. 


Example  5.3.  In  this  example,  the  body  forces  are  nonzero.  The  parameters 
for  this  example  are  E+  =  150,  E~  =  10,  v+  =  .2  and  v~  =  .24.  Figure  (7)  and 
Figure  (8)  are  the  mesh  and  contour  plots  of  the  displacements  obtained  from 
the  finite  element  method  developed  in  this  paper. 

From  these  examples,  we  see  that  the  numerical  results  obtained  from  the 
immersed  finite  element  method  agree  with  the  results  by  the  finite  difference 
method  in  [16,  17].  But  the  system  of  linear  equations  obtained  from  the  fi- 
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Figure  8:  The  mesh  and  contour  plots  of  the  computed  solution  v  of  Exam¬ 
ple  5.2  obtained  by  the  immersed  finite  element  method  developed  in  this  paper. 


Figure  9:  The  mesh  and  contour  plots  of  the  computed  solution  u  of  Exam¬ 
ple  5.3  obtained  by  the  immersed  finite  element  method  developed  in  this  paper. 


Figure  10:  The  mesh  and  contour  plots  of  the  computed  solution  v  of  Ex¬ 
ample  5.3  obtained  by  the  immersed  finite  element  method  developed  in  this 
paper. 
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nite  element  method  has  better  condition  than  that  obtained  from  the  finite 
difference  method. 


Conclusion  and  acknowledgments 

In  this  paper,  we  have  developed  the  immersed  finite  element  method  based  on 
a  uniform  Cartesian  mesh  for  elasticity  equations  with  discontinuous  physical 
parameters  across  an  interface.  The  interface  does  not  have  to  be  aligned  with 
the  mesh.  The  main  idea  is  to  modify  the  basis  function  over  those  triangles 
in  which  the  interface  cuts  through  so  that  the  natural  interface  conditions  are 
satisfied.  The  standard  linear  basis  functions  are  used  for  other  triangles.  A 
level  set  function  whose  zero  level  set  represent  the  interface  is  used.  Numerical 
examples  presented  here  show  second  order  convergence  for  the  new  method. 

The  first  author  was  partially  supported  by  NSF  grants  DMS-0201094  and 
DMS0412654,  and  by  an  ARO  grant  39676-MA. 
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